FVCOM_R_Subsetting

Using Bounding Polygons to Trim FVCOM Mesh Data

Published

March 25, 2024

Accessing FVCOM in R

This approach will lean on the work of Ben Tupper at Bigelow Labs, and their package fvcom. Many thanks to Ben and their team for their work on this.

Read more about their package here: https://github.com/BigelowLab/fvcom

Accessing NECOFS GOM4

As explained on the repo documentation for Ben’s {fvcom} package, Data are served via OpeNDAP on a THREDDS Server. NECOFS has its own THREDDS Directory which can be used to browse what is available.

Working with the Nodes & Elements

We will use this Jan, 2019 data to get information about the Nodes we wish to gather data from.

[1] "2019-01-01 00:00:00 UTC" "2019-01-01 01:01:52 UTC"
[3] "2019-01-01 01:58:07 UTC" "2019-01-01 03:00:00 UTC"
[5] "2019-01-01 04:01:52 UTC" "2019-01-01 04:58:07 UTC"

Working with the Mesh

Finding the Nodes We Care About:

We can do ourselves a favor by identifying which nodes we’re interested in before pulling variables. There are a couple different approaches to try:

  1. Use the nodes as points, use st_intersect to identify which are within the areas we are studying.
    2. Use the mesh itself, to preserve the triangular geometries. Then use the {fvcom} functions to pull data

Areas We are Interested in

  1. VTS Survey, Nearshore Lobster Habitat
  2. GOM & GB Ecological Production Units
  3. Northeast Shelf

Mesh Containment Within an Area

Will test this first with the offshore areas. But plan is to take the mesh and check for overlap/containment within the above polygons

Coordinate Reference System:
  User input: EPSG:6318 
  wkt:
GEOGCRS["NAD83(2011)",
    DATUM["NAD83 (National Spatial Reference System 2011)",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Puerto Rico - onshore and offshore. United States (USA) onshore and offshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore and offshore."],
        BBOX[14.92,167.65,74.71,-63.88]],
    ID["EPSG",6318]]

Use New Mesh to Extract Vars:

Extracting Vector Elements

Interpolate Currents to Regular Grid

Next Steps - GOM3

To go back further than 2016 we need to use the gom3 hindcast records. These have a different file location on the THREDDS directory, and we’ll need to confirm the grid structure and variables are consistent across these different model runs.

Verify that the IDS for nodes and elements are the same between NECOFS and FVCOMGom3.

If we want to go back further in time then NECOFS won’t yet be online and we’ll need to use data from the FVCOM GOM3 hindcast. Best-case scenario these mesh id’s match so we can append the time-periods together

Can we verify they match? Good.

What about the time frequency, what is the pace of the hindcast data records? Looks hourly, cool thats probably too much.

[1] 31.26466

Save Node/Element ID’s for regions of interest

Now that we know that the node and element ID’s are consistent across GOM3 & NECOFS, we can use the approach above to save the mesh information for areas within the regions we are looking at:

These can be used if-needed to limit the data being transferred over the OPeNDAP connection and save us data transfer/processing time. Would also be useful to have just as a means to match points to a standalone mesh geometry.

Coordinate Reference System:
  User input: EPSG:6318 
  wkt:
GEOGCRS["NAD83(2011)",
    DATUM["NAD83 (National Spatial Reference System 2011)",
        ELLIPSOID["GRS 1980",6378137,298.257222101,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["Puerto Rico - onshore and offshore. United States (USA) onshore and offshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore and offshore."],
        BBOX[14.92,167.65,74.71,-63.88]],
    ID["EPSG",6318]]

Download data for an extended period of time.

We don’t need hourly, or even daily records. But we will need to make sure we’re aggregating in time correctly and not just grabbing snapshots.

So where we sit now. We know the exact ID’s for the nodes and centroids that make complete mesh elements that fall within our area of interest. We aslo know the variables we are interested in, and possibly the time period we want to cover and the temporal resolution we’re interested in (months).

$names
 [1] "name"          "len"           "unlim"         "group_index"  
 [5] "group_id"      "id"            "dimvarid"      "units"        
 [9] "vals"          "create_dimvar"

$class
[1] "ncdim4"
$names
 [1] "name"          "len"           "unlim"         "group_index"  
 [5] "group_id"      "id"            "dimvarid"      "units"        
 [9] "vals"          "create_dimvar"

$class
[1] "ncdim4"
$names
 [1] "nprocs"        "partition"     "x"             "y"            
 [5] "lon"           "lat"           "xc"            "yc"           
 [9] "lonc"          "latc"          "h"             "nv"           
[13] "nbe"           "ntsn"          "nbsn"          "ntve"         
[17] "nbve"          "a1u"           "a2u"           "aw0"          
[21] "awx"           "awy"           "art2"          "art1"         
[25] "iint"          "Itime"         "Itime2"        "Times"        
[29] "zeta"          "file_date"     "u"             "v"            
[33] "omega"         "ww"            "ua"            "va"           
[37] "temp"          "salinity"      "km"            "kh"           
[41] "kq"            "q2"            "q2l"           "l"            
[45] "short_wave"    "net_heat_flux" "uwind_stress"  "vwind_stress" 
[49] "fvcom_mesh"   
$hasatt
[1] FALSE

$value
[1] 0
$hasatt
[1] FALSE

$value
[1] 0
[1] 48451    46
[1] -1.000001  0.000000
[1] 99137
start count 
    1    99 
[1] 6